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Gravitational waveforms which describe the inspiral, merger and ringdown of coalescing binaries 
are usually constructed by synthesising information from perturbative descriptions, in particular 
post-Newtonian theory and black-hole perturbation theory, with numerical solutions of the full 
Einstein equations. In this paper we discuss the “glueing” of numerical and post-Newtonian wave¬ 
forms to produce hybrid waveforms which include subdominant spherical harmonics (“higher order 
modes”), and focus in particular on the process of consistently aligning the waveforms, which re¬ 
quires a comparison of both descriptions and a discussion of their imprecisions. We restrict to 
the non-precessing case, and illustrate the process using numerical waveforms of up to mass ratio 
q = 18 produced with the BAM code, and publicly available waveforms from the SXS catalogue. 

The results also suggest new ways of analysing finite radius errors in numerical simulations. 
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I. INTRODUCTION 

Coalescing compact binaries are among the most 
promising candidates for a first direct observation of the 
gravitational waves (GW) predicted by Einstein’s Gen¬ 
eral Relativity. Within the next few years a new genera¬ 
tion of detectors composed by Advanced LIGO [TH3] , Ad¬ 
vanced VIRGO 0 and KAGRA [5], are expected to come 
online with design sensitivities ten times larger than the 
previous generation (LIGO and VIRGO). The detection 
and accurate identification of compact binary coalescence 
(CBC) events relies on the accuracy with which the ex¬ 
pected signals are modeled. The early inspiral stage of 
a CBC is well modeled by means of the post-Newtonian 
(PN) approximation [Bj, which consists of a weak-field 
slow-motion expansion of General Relativity in powers of 
GM/rc 2 , or equivalently u 2 /c 2 . However, as the binary 
tightens, the PN approximation becomes less accurate. 
In order to model the late inspiral one needs to solve the 
full Einstein equations, which in turn requires the use 
of numerical methods. An overview of the capabilities 
and techniques of numerical relativity (NR) as applied 
to black hole coalescence is given in BIB]. When the 
final state is a black hole, the perturbed Kerr black hole 
resulting from the merger will settle down emitting expo¬ 
nentially damped gravitational waves known as quasinor¬ 
mal modes, whose complex frequencies can be computed 
analytically in linear perturbation theory [9], while the 
amplitudes have to be determined from NR. 

In order to synthesize a GW signal that represents 
all of these stages, typically referred to as an Inspiral- 
Merger-Ringdown (IMR) waveform, one often needs to 
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align and appropriately glue together the PN and NR de¬ 
scriptions. The construction of such “hybrid waveforms” 
for the dominant quadrupole modes (corresponding to 
the i = 2, \ m\ = 2 spin-weighted spherical harmonic), 
which is what CBC searches have employed to date, has 
been discussed extensively in the literature, performing 
the construction both in the time and Fourier domains, 
see e.g. mm- A publicly available catalogue of such 
hybrid waveforms is described in [7j, their use for cali¬ 
brating search codes by injecting them into noise in m- 

The present paper focuses on complete waveforms 
that also include spherical harmonic modes other than 
£ = 2, \m\ = 2, or “higher order modes” (HOM). Our 
construction will be performed in the time domain (for 
work in the frequency domain see mm)- For simplic¬ 
ity we will restrict ourselves to non-precessing binaries, 
i.e. when the orbital plane is preserved. Hybrid wave¬ 
forms are essentially constructed in two steps. First we 
need to align the PN and NR waveforms, taking care of 
ambiguities of waveforms that describe equivalent physi¬ 
cal systems. Then we need to “blend” both waveforms to¬ 
gether. It is in particular the process of alignment where 
the presence of more than one spherical harmonic adds 
subtlety. The aims of this work are to provide a construc¬ 
tion algorithm for waveforms used in injection studies or 
waveform models including HOM, as well as to evalu¬ 
ate the fidelity of the input and resulting hybrid wave¬ 
forms, e.g. for constructing waveform catalogues similar 
to |7], which include higher modes. The type of wave¬ 
forms constructed here have been already used in [14] for 
the purpose of testing the effect of higher order modes in 
non-spinning CBC searches. We here present a detailed 
procedure for such a construction and analyze the main 
new sources of error that might appear. Our findings 
regarding the £ = 2, |m| = 2 modes are consistent with 
previous work, where the influence of several ingredients 
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of the construction (e.g. the choice of the hybridization 
point [interval], and the shape of the blending functions) 
has been studied extensively, see e.g. imiiaiisi. 

Our numerical waveforms have been taken from the 
publicly available SXS catalogue jT7] (computed by the 
SpEC code [MM]), and from a set of waveforms that 
have recently been constructed with the BAM code 
[30 EH- Our examples will focus on mass ratios 8 and 
18, where contributions from higher modes are signifi¬ 
cantly stronger than for roughly equal masses, but where 
it is computationally much more expensive to extend NR 
calculations to low frequencies, where PN is reliable. 

The rest of this paper is organized as follows. In 
Sec. [II] we introduce the basic definitions regarding wave¬ 
forms and their ambiguities; the basic definitions regard¬ 
ing data analysis which we will need are introduced in 
Sec. |III| In Section [IV[ we discuss the general principles 
of the construction of hybrid waveforms. In particular 
we discuss the three degrees of freedom required to align 
both waveforms. Section m is devoted to the review of 
the construction of a single-mode hybrid and the corre¬ 
sponding sources of error are identified and discussed. In 
Section lVD we describe the method we followed for coher¬ 
ently constructing our hybrids with higher modes and in 
Section [Vll] we study the residual disagreements between 
PN and NR in the hybridization region after alignment 
and identify their origin. Finally, in Sec. |VIlI| we evaluate 
the influence of NR extraction radius (and extrapolation) 
on the hybrid in terms of waveform matches. 

II. WAVEFORM DEFINITIONS AND 
AMBIGUITIES 

Our goal is to assemble a hybrid waveform from two in¬ 
dependently constructed pieces (the early-time part typi¬ 
cally computed by post-Newtonian methods and the late¬ 
time part computed by solving the Einstein equations nu¬ 
merically), or more generally to compare any two wave¬ 
forms (such as the results of two numerical relativity cal¬ 
culations) . As a first step we then need to understand the 
different conventions and possible ambiguities that went 
into the definition of both pieces. We will start defining 
our waveforms in terms of the Newman-Penrose scalar 
d' 4 , which is the waveform quantity directly computed 
in many numerical relativity codes, and afterwards focus 
on the strain h, which is the quantity directly relevant to 
the data analysis of current ground-based gravitational 
wave detectors. 'P 4 is computed by contracting the Weyl 
tensor C atl pv with the appropriate elements of a suitable 
null tetrad to m , to p , (see [52 for a detailed de¬ 
scription of the formalism). As an example, the BAM 
code uses the definition 

T 4 = -C ai ip„n^n v m a m 13 , (1) 

where i and n are ingoing and outgoing null vectors and 
-£-n= 1 = m-fh. The precise choices made in this code 
can be found in Sec. Ill of [501). 


The definition of 4^ carries with it several ambiguities, 
starting with the overall sign convention for the Riemann 
and Weyl tensors (including metric signature). As exam¬ 
ples, in the BAM code [50] the conventions from Misner, 
Thorne and Wheeler [55] are used, and the opposite sign 
is used in the SpEC code (see e.g. the comment above 
Eq. (2.100) of [M]). In addition, the overall sign in (TJ is 
a convention that may change between different authors. 

Furthermore, there is freedom in the choice of the 
tetrad. While n^ will coincide between different codes 
in the limit r —> oo, there is no canonical choice [55] 
for the complex null vector m^ which can be rotated by 
some angle er (m M —> e l<T in^), leading to a redefinition 
\l/ 4 —> e -2 * <T 'I' 4 . The different choices in the definition of 
T 4 thus amount to an ambiguity \l/ 4 —> e*^°'I' 4 , which in 
physical terms is simply the freedom in defining the two 
gravitational wave polarizations. 

The two real polarizations h + and h x of a gravitational 
wave can be conveniently represented as a complex strain 

/i(t, 9, ip ; S) = h + (t, 9 , p- S) - ih x (t, 9 , ip ; S), (2) 

where t is an inertial coordinate at null infinity, 9 is cho¬ 
sen as the angle between the line-of-sight from the detec¬ 
tor to the source and the total angular momentum of the 
binary (which we choose as our z-axis), ip is an azimuth 
angle in the source frame, and the intrinsic parameters 
of the source are collectively denoted as 5. This quan¬ 
tity can be obtained from Tij by applying a double time 
integration (see [35] for a discussion of the issues arising 
in this procedure), or directly from projecting the metric 
perturbation onto some some orthonormal polarization 
triad as is usually done in the PN context. Different 
choices of triad will again lead to a redefinition of the 
type h —> e l ^°h (see for instance Eq. (2.6) of [57]). 

It is convenient to decompose the strain into spin —2 
weighted spherical harmonic modes as 

OO i 

h(t, e, ip-, 3) = ( 3 ) 

1—2 m=-£ 

where Yr'f (0. <p) are the spin -2 weighted spherical har¬ 
monic basis functions. 

We restrict to the non-precessing case, here the in¬ 
trinsic parameters are then the total mass, the mass ra¬ 
tio and the two (dimensionless) spin projections, S = 
{M, q, XI 1 X 2 } onto the angular momentum of the sys¬ 
tem. In this case the geometry is symmetric with respect 
to the orbital plane, which is preserved in time. This 
equatorial symmetry implies 

h{t , 6, tp\ S) = h * (i , 7T - 9, ip-, S) (4) 

(where a * denotes complex conjugation) provided that 
the polarizations are defined using some appropriate 
choice for the projection triad/tetrad, which is usually 
the case in the literature. For the individual modes this 
translates into 

he,m(t, 2) = (-1 ) f /i|_ m (t,S). 


(5) 
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Therefore, we just need to focus on the m > 0 modes, ex¬ 
cept when reconstructing the whole waveforms. Finally, 
it is convenient to decompose each mode into a real am¬ 
plitude and phase as 

S) = A e>m (t, (6) 


In the following we will omit the dependence on S in 
order to simplify notation and write h(t, 9, p). 

We note that during inspiral the phase of the (£, m) 
mode approximately follows the rule « m0 OT \,(t), 

where (p OI b is the orbital phase, however this approximate 
relation has to break down eventually during the merger, 
as it is violated during the ringdown, as one can check by 
comparing the quasi-normal frequencies of the different 
modes. We will return to this issue in Sec. IVII Bl when 
comparing PN and NR phases of individual modes. 

The strain h D seen by a detector located in the direc¬ 
tion (9, p) of the source sky also depends on the luminos¬ 
ity distance to the source, and the orientation of the 
detector with respect to the source, which we parametrize 
using three angles (8,0,0). Here (0,0) are the spherical 
coordinates of the source in the detector sky, and 0 is 
a polarization angle. This dependence is encoded in the 
antenna patterns F + and F x of the detector as 


D _ F+(0, ip, 0)h+(t, 9, p) + F x ( 6 , ip, 0)h x (t, 8 , p) 

d L 

(7) 

where 

1 -f- cos^ 9 — 

F + = ---cos 20 cos 20 — cos 9 sin 20 sin 20, 

1 -f- cos^ 9 — 

F x = ---cos 20 sin 20 + cos 9 sin 20 cos 20. 

It is well known that this can be rewritten as 


h D = ^(0,p, 4 ) [ cos K } l+ [t, 9,p) -f- sin k h x (t, 0, p)] , 

«L 

(8) 

where n acts as an effective polarization angle and F/d^ 
is a simple overall amplitude factor. 

We can now list the possible ambiguities in the defini¬ 
tion of the waveform and its spherical harmonic modes for 
two waveforms A and B, computed with different meth¬ 
ods and conventions. We use the superscripts A and 
B to refer to quantities derived from these waveforms. 
For aligned-spin binaries we assume that computations 
A and B preserve the manifest equatorial symmetry of 
the problem, in particular that the z-axis of the coor¬ 
dinate system we use to define our spherical harmonic 
mode decomposition is parallel to the angular momen¬ 
tum of the system. The remaining conventions to choose 
are the origin of the azimuthal angle p of the spherical 
coordinates, a polarisation angle 0o, and the origin of 
the time coordinate. Neglecting for the moment issues 
related to the accuracy of computations A and B, the 
strains h A and h B computed by implementations A and 
B are related by 

h A (t, 9, p) = e l ' l ’ 0 h B (t + r,9,p + p 0 ), (9) 


where 0o and po are two angles that encode the differ¬ 
ent choices in conventions. Of course, the same relation 
applies to 'F 4 . As a result, the /q. m modes are related by 

htm(t)=e i ^ +m ^hf tm (t + r) ( 10 ) 

Usually, conventions are chosen such that Eq. ([ 5 ]) holds 
for the individual modes. This implies that 0o £ {0, it} 
and thus 


h A m(t)=(-lY°e im ^h B m (t+T) ( 11 ) 

with Ko £ {0,1}. In the case where one only considers 
the dominant (2,2) mode, equations (101 and (11) can 
be rewritten as h A 2 (t) = e l ^°h B 2 (t + r) i.e. the whole 
angular freedom amounts to a global phase shift. In the 
multi-mode case, comparing waveforms without ensuring 
a consistent choice of 0q will lead to incorrect results. 

In order to compare or hybridize two waveforms, we 
thus need to align them to resolve the above ambiguities, 
i.e. either keep track of the differences in conventions 
(which is typically not possible) or infer these from the 
waveforms themselves. We describe in detail our pro¬ 
cedure to do so in Sec|VI[ which involves defining some 
notion of distance between the two waveforms and min¬ 
imizing it over the parameters (0o,po , r). Note that 0o 
only depends on the difference in the definition of the po¬ 
larizations between methods A and B and therefore needs 
to be determined only once whereas r and po will differ 
for each pair of waveforms. Inaccuracies in the waveforms 
will in general lead to residual discontinuities between the 
spherical harmonic modes even after alignment and these 
are studied in detail in Sec. m 


III. MATCHED FILTER ANALYSIS 


Gravitational wave signals buried in the noise n(t) of a 
detector can be extracted using the technique of matched 
filtering, which is the optimal filtering to extract signals 
of known shapes buried in stationary Gaussian noise. 
Much of the complications of gravitational wave data 
analysis with actual data from detectors arise from non- 
stationarity and in particular non-Gaussian noise contri¬ 
butions, which we will not consider in this paper. There 
is a variety of effects which may cause a calculated wave¬ 
form to differ from the physically correct one. The rele¬ 
vance of such differences will depend on the application 
of the waveforms. It is well known that inaccuracies in 
the template will in general degrade the performance of 
the matched filter, and signals may be missed or their 
parameters (such as masses and spins) incorrectly identi¬ 
fied. Our analysis of the match of waveforms will follow 
the lines of the matched filtering procedure for detecting 
gravitational waves. We now briefly review some of the 
data analysis concepts which we will use. 

For time domain real waveforms h\(t) and h 2 (t) one 
defines the inner product (“Wiener scalar product”) 


</n | h 2 ) = m 


hijmui 

Snif) 


0 


( 12 ) 
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where S n (f) is the one-sided power spectral density of 
the noise n(t), and tildes denote the Fourier transform of 
the respective time series. In defining this inner product, 
we use the fact that the Fourier transform of the real 
gravitational wave strain h(t) satisfies h(f) = 
thus allowing us to define the inner product as an integral 
over positive frequencies only. In practice, this integral 
will be performed in a frequency band determined by 
the detector bandwidth, the lower cutoff frequency fo 
being given by the seismic wall of the detector noise. In 
this paper we consider the predicted 2015-early science 
run noise curve for Advanced LIGO with a 30 Hz lower 
frequency cutoff [3B], as used in El. as well as its optimal 
sensitivity given by the “zero-detuned high-power” noise 
curve [39] with a 10 Hz lower frequency cutoff. 

The overlap between two waveforms is defined as their 
normalized inner product 


oi^m 


(hi | h 2 ) 

\J(hi | hi)(h2 | / 12 ) 


(13) 


and the match as the maximized overlap over a chosen 
set of extrinsic parameters {A*} of the waveform, 


M[hiM] 


maxOf/ii, 

{A,} 


(14) 


The extrinsic parameters {A,} may be, for instance, the 
time of arrival and the coalescence phase of the binary 
{i 0 ,^o}- As an example, the maximization over to is to 
be understood as the maximum overlap between hi and 
/i 2 obtained by shifting in time one of the templates, i.e., 
h(t) —t h(t + to). In the case of waveforms containing 
only a single pair of (£, |m|) modes, e.g. the dominant 
( 2 , 2 ) and ( 2 , — 2 ) modes, the maximization over all the 
extrinsic parameters can be performed analytically as 


M[hi,h 2 ]= max 0 [hi,h 2 ] = 


= max 

to 


r°° hi U)h% (/) „-i 27 T ft n 
Jfo S„(f) 6 aJ 


\J(hi | hi)(h2 | ^ 2 ) 


(15) 


whereas in general the maximization over 9, ip, and if) 
needs to be done numerically. We will however only con¬ 
sider maximisation over to, <po hr this paper. 


IV. NR AND PN INPUT WAVEFORMS 

1. Post Newtonian expansions 

PN expansions compute an approximate solution of 
Einstein’s equations up to a certain order in the expan¬ 
sion parameter. Since only a finite (small) number of 
expansion terms are known, it is not possible to perform 
a strict convergence test to estimate the truncation error. 
In addition, the approximation breaks down at merger or 
shortly before. Even at a given PN order for the energy 
and the flux, different treatments in the derivation of the 


orbital phase from the balance equation give rise to a va¬ 
riety of “PN approximants”, such as the Taylor approxi- 
mants [H ED, which are commonly used in gravitational 
wave data analysis. 

The main consequence of the PN truncation error is 
a phase evolution which progressively deviates from the 
correct one as the binary evolves. This secular trend 
translates into the key source of error for the estimation 
of the time-shift r between PN and NR, as shown in 
Fig|2] where the secular trend is also shown in compar¬ 
ison with oscillations originating in residual eccentricity 
of the NR data. Since secular phasing errors in PN grow 
with frequency, it is desirable to hybridize at low fre¬ 
quencies, or equivalently with very long NR waveforms, 
to minimize such errors. Longer NR waveforms are how¬ 
ever significantly more expensive computationally. 

In our study we use the Taylor T1 and T4 approxi¬ 
mants including 3.5 PN non-spinning [5] and spin-orbit 
m and 2PN spin-spin m phase corrections, which we 
will just denote as T1 and T4 for brevity. We use 3PN 
non-spinning amplitude corrections for the higher modes 
[44] and 3.5PN for the 22 mode [45]. The spin corrections 
to the amplitudes are known up to 2PN US). 


2. Numerical Relativity 


In NR, the expansion parameter of the PN approx¬ 
imation is essentially replaced by an expansion in the 
resolution of the numerical grid, and at least in principle 
it is possible to provide error estimates through a con¬ 
vergence test. Unless the GW signal is calculated at null 
infinity, e.g. employing the Cauchy characteristic method 
G23 0Z!, convergence needs to be checked not only with 
respect to grid resolution but also extraction radius. It 
is then possible to either extrapolate to infinity from a 
series of finite radii (e.g. for the SpEC simulations mb 
ED El EH]), or to directly use results from a single finite 
radius. Furthermore, systematic errors arise due to im¬ 
perfections of the initial data and initial orbital param¬ 
eters, in particular residual eccentricity, which generates 
oscillations in amplitude and phase (see e.g. the discus¬ 
sion in mi). Unphysical radiation content of the initial 
data manifests itself in a small GW burst at early times, 
which is usually referred to as “junk radiation”. 

For mass ratio q = 18, new NR simulations have been 
performed with the BAM code mi ed, and are sum¬ 
marized in Sec. III. A of mi- GW strain is computed 
from T 4 using the fixed-frequency-integration algorithm 
described in [35]. For a recent comparative discussion 
of current numerical relativity codes for black hole bina¬ 
ries, including SpEC and BAM see IDE]- Fig .[ljshows the 
amplitude of the most important modes for the highest 
resolution BAM q = 18 waveform. 
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BAM q=18,*=0 



t/M 


which implies uj B (t) = u) A (t + r). Determining r and po 
is then trivial: one chooses any frequency ujq inside the 
range [w B (0), u) A (t A )] and obtains 


r = t A {u 0 ) - t B (uj 0 ), e lcp ° 


h B (t B (uj 0 )) 
h A (t A (L0 0 ))' 


(17) 


In this idealized case the time alignment and angle p>o do 
not actually depend on the frequency coq and no blending 
is required, both functions perfectly overlapping before 
and after the matching point. 


B. Realistic case 


FIG. 1. Logarithmic plot of the mode amplitudes for the non- 
spinning (?=18 configuration. The modes with (£, m = l) (in 
black, with l = 3,4, 5 from top to bottom) and {l,m = £ — 1) 
(in dashed blue, with i = 2, 3,4 from top to bottom) all have 
a peak amplitude smaller than that of the (2,2) mode (on top) 
by a factor between ~ 3 and 20. We show the clean part of 
the waveform after initial transients due to junk radiation. 


V. SINGLE MODE HYBRIDS 
A. Idealized case 

In order to illustrate some key points of the hybrid con¬ 
struction, we first consider a single mode and assume that 
we have at our disposal two infinitely accurate general- 
relativity computations of some spherical harmonic mode 
of the strain, h A (t ) and h B (t ), that overlap over some 
portion of the evolution of the binary, i.e. that there is 
an interval where they satisfy ( [To] ) . Let us define for con¬ 
venience the amplitude A(t) and the phase <j>(t) of wave¬ 
form X (with X= A or B) as h x (i) = A x (t)e l( ^ W which 
we assume to be defined over some interval [0, t x ], as well 
as the frequency u x (t) = d<p x /dt which is a monotonic 
function of t in the case of binaries on circular orbits (we 
will discuss below the problems introduced by the oscil¬ 
lations due to some residual eccentricity in the NR wave¬ 
form; however, provided that the eccentricity is small 
enough, this remains true). We can therefore define the 
inverse function t x (u>), which satisfies t x (co x (t)) = t. 

Our idealized infinitely accurate waveforms will satisfy 

h B (t) = e iVo h A (t + T) (16) 

for some r and ipo, and t in the interval [t A (uj B (0)),t A ], 


In practice both computations are affected by errors, 
and (16) is never exactly satisfied over any interval. One 


rather has to find the best parameters r and ipo so that 
(161 is the closest to being satisfied in some sense and 


over some matching window. We thus now have to make 
some particular choices in our hybrid construction. We 
parameterize our window by the initial time to or initial 
frequency ujq (defined as uj PN (t 0 ) = oj o) and the length 
of the window in the time domain At, i.e. our window is 
[to, to + At]. In order to avoid the influence of amplitude 
errors, we only take into account phase information when 
aligning the waveforms in time. We adopt the quantity 

rto+At 

A(r; t 0 , At) = / (u, PN (t) - co NR (t r)) dt, 

Jt 0 

(18) 

which has the advantage of not depending on po■ Other 
authors have replaced to by (/> or h [15] in the integrand 
(which then depends on <po) and tested that their hybrid 
was not affected much by this choice. 

The appropriate time shift r for a given choice of 
window [to, to + At] is then obtained by minimizing 
A(r;to,At). Once this is done, the optimal phase shift 
ipo has to be determined. Simple choices are to align the 
phases at a fixed time, e.g. the beginning of the window, 
ip 0 = 4> NR {t o — r) — (j) PN {to), or to pick the phase shift 
that minimizes J* 0+At — r) — <fi PN (t) + v?o)~ dt. 

We have checked tirat the resulting hybrid depends very 
weakly on this particular choice. This is due to the fact 
that the phase has one additional integration with respect 
to the frequency, so it contains less oscillations. 

Once r and po have been determined, both waveforms 
are combined into a piecewise definition 


h{t) 


e iv °h PN {t + r) 

w-{t)e^°h PN (t + t) + w+(t)h NR (t) 
h NR (t) 


t < to — t 

if to — t < t < _ T + At 

if to — t + At < t 


(19) 
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where, with the notation (t) for blending func¬ 

tions that monotonically go from 0 to 1 (or from 1 
to 0) in the interval [ti, < 2 ]? we have defined w ± (t) = 
w t 0 _ T t 0 _ r +At](^)’ i- e - we P er f° rm the blending over the 
same interval we used to determine r and ipo- Here again, 
different authors have made different choices for the ex¬ 
act shape of these functions. For instance, m considers 
cosine functions, while we here use linear ones. 

Note that multiplying the PN part by e l<Po is a redef¬ 
inition of the conventions for the PN waveform (change 
of the orbital phase): now the early part of the hybrid 
does not exactly reduce to the original PN waveform. In 
the single mode case this is trivial (and multiplying the 
NR part by e~ ltpo would have been equivalent). As sin¬ 
gle mode matches are always optimized over coalescence 
phase, this redefinition of conventions will never have any 
practical consequence. As we will see, in the multimode 
case, things get more involved. 



- A<A=ti/2 

- A0=27T 

- A0-47T 

- A0=5;r 

- A0=7;r 





FIG. 2. Top: relative time shift (r — tq)/M as a function 
of Mujq for q = 3 non-spinning SXS NR data hybridized to 
T1 for various lengths A (f> of the matching window. The ref¬ 
erence time shift to is the one obtained for A (f> = 7n and 
Moj 0 = 0.043. Note how for larger A(j> the oscillations in the 
estimation of r/M are smaller. Bottom: same for several SXS 
data sets hybridized to T1 and T4 using A <f> = 7n (smaller 
values of q on top) and to as above. 

In order to quantify the error in the time alignment of 
the waveforms, it is instructive to look at the time shift 
r as a function of matching frequency uj 0 and window 
size At (note that the absolute value of r for some ujq is 
meaningless, what matters is its variation). Figj2] illus¬ 
trates how our best choice for r varies with our choices of 
window length At, and how the secular trend depends on 


the choice of PN approximant. In the present case, oscil¬ 
lations in the NR waveform are caused in particular by 
residual eccentricity, which manifests itself at frequencies 
of the order of the orbital frequency (close to half the fre¬ 
quency of the (2,2) mode). Indeed, for At significantly 
larger than the orbital period, we see that most of the 
oscillations in the NR data average out, and for a At cor¬ 
responding to at least Ac/) = cj) PN (t 0 +At)—(j) PN (t 0 ) ~ 57 t 
oscillations are smaller than the secular trend due to the 
phase evolution not being accurate. Unless specified oth¬ 
erwise, we therefore choose At such that A cf> = 7 tt, i.e. 
3.5 GW cycles. Another possibility, proposed in [511] is to 
force the window extremities to lie at some maximum of 
the modulation to ensure the cancellation of the effects 
due to the modulation over the window. 

In order to minimize alignment errors due to the sec¬ 
ular dephasing between PN and NR, the interval over 
which one aligns the waveforms should be chosen as early 
as possible since the accuracy of the PN perturbative 
treatment degrades as the frequency increases, but not as 
early as to be affected by junk radiation or other early- 
time transient errors. In addition, a comparison of PN 
approximants as in Fig|2] can be used to choose a PN 
waveform with smaller error in the matching region. 



M[M 0 ] 


FIG. 3. Effect of At on the final hybrid (2,2) mode. We 
show the match A4 [h TQ , h& T \ optimized over time and phase 
between hybrids built from SXS q = 3 non-spinning data 
hybridized to T1 at Muj 0 = 0.043 (solid) and 0.073 (dashed). 
An artificial time-shift to + At is applied when constructing 
Ha t , while h TQ has been built using the optimal time-shift to 
between PN and NR. See main text for details. 

Now we address the question of how much this im¬ 
pacts the waveform in terms of quantities useful for data- 
analysis. Fig. [3] shows the match between q = 3 non¬ 
spinning hybrid waveforms constructed using artificial 
time shifts r = t 0 + At with a reference waveform for 
which t = To- We use the Zero-Detuned High Power 
noise curve of Advanced LIGO m with a lower frequency 
cutoff of 10Hz in order to facilitate comparison with HD- 
Regardless of the intrinsic parameters of the system and 
the hybridization frequency, the mismatch increases with 
At. When the hybridization region is in band, the match 
decays by a few 10 -3 for a value of At of a few M , which 
is consistent with the results obtained by [15] . Note that 
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At has a larger effect on the match for larger frequency 
wo, as expected from the fact that At is then a larger 
fraction of the period. 


VI. MULTI-MODE HYBRIDS 

In this section, we describe our procedure to construct 
hybrid waveforms with higher modes and define quan¬ 
tities that will be used in the error analysis of Sec |VII| 
Higher modes become increasingly relevant for binaries 
with large mass ratio, for this reason we illustrate our 
procedure using a waveform produced recently with the 
BAM code for a non-spinning binary with q = 18, which 
we have presented in Sec. m 


A. Step 1: Determination of r, i/>o,y>o 


and the length of the time-window over which the wave¬ 
forms are aligned in time as At. Considerations on how 
to choose these two parameters have been described in 
the previous section. We can then first determine r by 
minimizing the same quantity as in the single mode case, 


rto+At 

A(r; to, At) = / (t) - w^ 2 fl (f - r))“ dt, 

Jt 0 

( 22 ) 

since only the frequencies enter in A(r;to,At) so that 
the determination of r decouples from that of the phase 
offsets (ipo, Vo)- Before moving on to the determination 
of (ip o, </?o), let us recall that given a code to generate PN 
waveforms and an NR code, (fo will depend on choices 
made to generate each individual waveform whereas ipo 
could in principle be computed once and for all by com¬ 
paring all the convention choices in both codes. In what 
follows, we assume that ipo can only take the values 0 or 
7 r as discussed in Sec. m Let us define 


As discussed in Section [TT| in the presence of higher 
modes one needs three parameters (r, ipo, ipo) to describe 
the possible differences in conventions between the PN 
and the NR waveforms. The best choice of such parame¬ 
ters will depend on the matching time or frequency, but 
it appears fruitful to not make different choices for dif¬ 
ferent directions in the source sky (it seems conceivable 
but not practical to do this ). Several strategies to infer 
these parameters from the waveforms can be explored. 
One important ingredient is how to weight the contribu¬ 
tion of different modes in determining (t, -ipo, Vo)- One 
natural choice, pursued in m, is to define a single set of 
(T,ipo,ipo) by minimizing the quantity 


A&.m = - T) ~ <^m(*0)- (23) 


Then ideally (i.e. assuming that 
ipo + 2</? 0 + A(p 2,2 = 0 mod 2n i.e. 


(10) holds), 


we have 


A(p 2 ,2 + i>o 

Vo = -x- 


mod 7T, 


(24) 


which gives 2 solutions for ipo in the interval [0, 2-7 t[ if ipo 
is previously known and 4 solutions if ipo is unknown but 
restricted to ipo £ {0,7r}: 


(ipo,Vo) 




2 > 


7 T mod 27T 


(25) 


/ diXXmt* - T )e^ +m ^ h™(t )|, (20) 

J e,m 

where the integral is performed over some window cor¬ 
responding to the hybridization region, and the contri¬ 
bution of each mode is naturally weighted with its am¬ 
plitude. Note that m does not restrict ipo to belong to 
the set {0,7r}, and the resulting modes do not in general 
follow the usual relation ([5]). 

In this paper, we take a different approach, constrain¬ 
ing the 3 degrees of freedom as much as possible using 
only the dominant (2,2) mode. The (2,2) mode of hybrids 
constructed this way will thus coincide with hybrids con¬ 
structed only for the (2,2) mode, and two hybrids con¬ 
structed with different sets of higher modes will exactly 
coincide on their common modes, which facilitates com¬ 
parisons and studies of the contribution of some specific 
mode. In practice, our procedure is as follows. Just as 
in the single mode case, we parametrize how early (or 
late) in the evolution we perform our hybridization using 
a “hybridization frequency” to o, which defines the “hy¬ 
bridization time” to through 

-^■(*o)=wo, (21) 


with k £ {0,1} and n' £ {0,1}. To lift this degener¬ 
acy, we need information from at least one of the higher 
modes, say (£*, to*), and we use the one with the largest 
amplitude, typically the (3,3) mode unless it is zero for 
symmetry reasons. If again both waveforms were in¬ 
finitely accurate, (101 would imply 


ipo + m*ipo + A <pt,m, = 0 mod 27T, 


(26) 


but in the presence of waveform errors this will not hold 
for any of our four solutions. However, we can choose 
the solution that is the closest to satisfying this equation, 
which is uniquely determined only in the case where ?n* is 
odd. Note that in the case where only even m modes are 
available, ipo needs in fact only to be determined modulo 
7 T since only the combination rrupo appears in the hybrid 
construction and the two solutions can be discriminated 
using any even higher mode. 

We find that there is a relative ipo shift of n between 
the BAM and SXS waveforms, and also between BAM 
and the conventions used in the PN context by Arun et 
al [37] and Blanchet [5j. 

In Fig [4] we show the solutions that were found using 
this procedure on our q = 18 case hybridized with three 
different PN approximants as a function of the hybridiza¬ 
tion frequency. Not surprisingly, the result corresponds 










— TaylorTI 

— TaylorT4 
- SEOBNRv2 


FIG. 4. Estimations for ipo as a function of Mujq for the 
case of q = 18 non-spinning BAM NR data hybridized to PN 
T4, T1 and SEOBNR.v2 (from top to bottom). For Tl, the 
estimation of tp o changes by ~ 7r over the frequency range 
shown here (corresponding to ~ 7 cycles in h, 2 p ). 


to the lower plot of Fig. [2] which show the dependence of 
the hybridisation time shift on frequency and PN approx- 
imant. In Fig [4] we see that both the standard Tl and 
T4 approximants exhibit large secular trends, indicating 
a large difference in the orbital phase (or more precisely 
the phase of the (2,2) modes) between the PN waveforms 
and the NR result over the frequency range shown here. 
In contrast (and remarkably), the SEOBNRv2 waveform 
m shows almost no secular trend even though the model 
was calibrated to NR simulations only up to q = 8. Eval¬ 
uating the secular trend of the PN approximant as com¬ 
pared to the NR waveform is an important part of the 
hybridization procedure. Exactly as in the single mode 
case, this secular trend translates into some “hybridiza¬ 
tion error” (for instance, the phase of the (2,2) mode for 
two hybrids built using Tl but with Mojq = 0.085 or 
Mu> o = 0.115 and aligned in the early inspiral will dif¬ 
fer at the peak by almost one gravitational wave cycle) 
but this error has nothing to do with the higher modes 
themselves and controlling it is not our main focus here. 
Instead, we will try to identify additional figures of merit 
for the hybrid that directly quantify the additional error 
due to the higher modes. 


B. Step 2: evaluate residual disagreement between 
PN and NR at the matching point 

We now investigate the residual phase and amplitude 
disagreements between PN and NR at the matching point 
and define appropriate quantities to describe this dis¬ 


{ e i(m V0 +1> 0 ) h PN(t + T } 

wi m (t)e i ^°+^h PN (t + t) + 
h» R (t) 

with t 0 £ + At). Fig. [fi] shows three of the 

resulting hybrid modes. 


agreement, while we postpone the analysis of the main 


source of this disagreement to Sec. VII 


In the idealized case, correcting for the differences in 
conventions using (ipo,ifo) is sufficient to ensure that all 
the phases are continuous between the PN and NR wave¬ 
forms at the matching point, i.e. that the quantities 


Q,m = A (t>^ m + 4>o + rrupo, 


(27) 


which are functions of the hybridization frequency ojq are 
all zero. In practice, this is not the case and we will use 
these quantities as measures of the residual phase dis¬ 
agreements. The values for our example q = 18 case are 
shown for the most important modes in Fig. [5] (left) and 
are typically of a few degrees for m = £ modes and 10 — 15 
degrees for m = £ — 1 modes. Apart from the values 
themselves, one important feature is the fact that unlike 
for ipo, these remain roughly constant over the range of 
hybridization frequencies explored here. This is a conse¬ 
quence of the fact that we are essentially measuring phase 
differences between the higher modes after aligning the 
(2,2) modes at coq via the term nup o in Eq. (271, which 
effectively absorbs the secular dephasing shown in Fig. [4] 
In other words, the really quantify the residual dif¬ 
ferences between PN and NR in the dephasings between 
the higher modes and the (2,2) mode, factoring out the 
error in tracking the orbital phase (or equivalently that of 
the (2,2) mode) of the system. Regarding the amplitude 
discrepancies, we can simply define the ratio 


Tim — 


l^m (*0 ~ T )l 

\hZ(to)\ 


(28) 


which is plotted in Fig. [5] (right). We find that modes 
with higher frequencies show larger amplitude disagree¬ 
ments. We will perform a detailed analysis of the phase 
and amplitude errors in Sec. |VII| 


C. Step 3: hybrid construction 


We can finally proceed to constructing the higher hy¬ 
brid modes as piecewise functions in a similar fashion 


as in Eq. (191, however we now allow different blend¬ 


ing functions for different modes to deal with noisier low 
amplitude modes. With the notation [ty 171 — T,t e 0 ’ m — 
t + At 1 ’™] for the interval used for each mode and 






(f) the associated blend¬ 


ing functions, we define 


t < t e 0 ' m — r 

wj m {t)h NR {t) 4 ’™ - t < t < t e 0 ’ m - t + Atf'™ (29) 

4’ m - r + A t e ’ m < t 


VII. HYBRID HIGHER MODES: PHASING 
AND AMPLITUDE ERRORS 


In the previous section, we introduced figures of merit 
to quantify the residual disagreements between the PN 
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FIG. 5. Left: Phasing errors ee, m (in degrees) for the q = 18 non-spinning BAM data hybridized with T1 as a function of Muj o. 
Note the lower values for l = m modes (black) as compared to i m (red). Right: Amplitude ratios rq m for the same hybrid 
construction. Note how modes with larger m, i.e. with larger frequencies, tend to show larger amplitude disagreements. 


and NR higher modes both for the phase (see Eq. {27}) 
and for the amplitude (see Eq. (281), after correcting for 
the differences in conventions by aligning the (2,2) modes. 
We devote the present section to identifying the main 
source of these disagreements among the errors affecting 
both computations and described in Sec. m In partic¬ 
ular, we find that the phasing errors ee m are dominated 
by the fact that the waves are extracted at finite radius 
in the NR simulations. Regarding the amplitude errors 
rim, we observe that extraction radius in NR simulations 
plays an important role (dominant for some modes), but 
that PN truncation errors are dominant for some other 
modes. For a detailed analysis of errors in higher modes 
for a simulation of non-spinning q = 4 numerical sim¬ 
ulations see m- One of the main conclusions there is 
that extrapolation of gravitational waveforms to infinite 
extraction radius is particularly important for subdomi¬ 
nant multipoles with t |m|. 


For this study, we have looked at a variety of physi¬ 
cal configurations (mass ratios, spins), and we focus on 
those configurations for which several NR simulations 
have been performed using different codes (and there¬ 
fore different numerical setups, gauge conditions and ini¬ 
tial data). We illustrate our results with the case of a 
q = 8 non-spinning binary simulated using the BAM 
code [SOI EH, and also available in the public SXS cata¬ 
logue m This is an interesting case with strong higher 
mode contributions due to the large mass ratio and where 
wave signals are available at different resolutions for both 
codes, as well as several extraction radii. Furthermore, 
the SXS data sets are also extrapolated to null infinity 
at different polynomial orders, see eu for a discussion of 
different methods, and m for a comparison with charac¬ 
teristic extraction results. Note that the SXS waveform 
is significantly longer, which allowed us to hybridize at a 
frequency as low as Muo = 0.043, when the BAM wave¬ 
form requires Mujq > 0.080. This means hybridizing 40 
and 9.2 gravitational wave cycles (in the 22 mode) before 
merger respectively. 


A. Amplitude errors 


We first focus on the residual amplitude discrepancies 
measured by the ratios ry„ ( defined in Eq. (28) and start 
by investigating the effect of finite radius extraction on 
these quantities. Fig.[T]shows the amplitude ratio for dif¬ 
ferent modes as a function of the hybridization frequency 
wo, and for different extraction radii for both BAM and 
SXS, including the SXS waveform extrapolated to null 
infinity with a polynomial of order N = 4. Note that all 
these curves are amplitude comparisons between NR and 
PN, but by taking the ratio of two curves, one obtains a 
direct comparison between two NR results. 


While the amplitude of the waves extracted at finite ra¬ 
dius typically differs from that of the extrapolated waves 
by of the order of one percent for the highest radii avail¬ 
able, significant errors arise when waves are extracted 
closer to the source, and errors are quite different for 
data sets computed with different codes. With the ex¬ 
ception of the (3,2)-mode, e.g. r = 100 BAM data are 
significantly closer to the extrapolated result than the 
corresponding SXS r = 100 curve. This is not surprising, 
since finite radius errors depend on gauge conditions, and 
the choice of lapse and shift are indeed different for both 
codes. We have at present no explanation why the BAM 
(3, 2)-mode shows anomalous behavior as compared to 
the other modes. We note that, for any given finite ex¬ 
traction radius, the error becomes larger at lower frequen¬ 
cies. This is the expected consequence of the fact that 
as the frequency of the waves increases (or equivalently 
as their wavelength decreases), the wavezone (defined by 
r )$> A) extends to smaller radii. 

Having highlighted the effect of extraction radius on 
the NR amplitude, we now come back to the comparison 
between NR and PN and restrict our attention to the ex¬ 
trapolated SXS amplitude on the NR side. For the (2,1), 
(2, 2) and (3,2) modes shown in Fig. [7] the ratio remains 
almost constant and differs from 1 by at most two percent 
over the whole range of frequencies considered. On the 
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FIG. 6. Amplitude (left) and real part (right) of non-spinning BAM q = 18 (2, 2), (2,1) and (3,3) modes hybridised with Tl, 
from top to bottom. We show PN (blue), NR (red) and hybrid (black) modes and focus on the hybridization region. The (2,1) 
mode is a typical example of good amplitude agreement and large Q, m , while the (3,3) mode is a typical example of small e^ m 
and poor amplitude agreement. Note the usage of different blending windows from the one used for the (2,2) mode. 


contrary, the ratio for the (3,3) mode shown in the lower 
right panel features a strong secular trend, significantly 
departing from 1 at high frequencies by a few 10%. We 
observed the same behavior in the other relevant higher 
modes (i.e. the (4,3), (4,4) and (5,5)) not plotted here. 
Some smaller (but still of the order of several percent) 
disagreements are visible at low frequencies. The agree¬ 
ment between the different NR curves (extrapolated SXS 
and the outermost extraction radius for BAM), at least to 
a much higher degree than the disagreement between PN 
and NR, and the fact that this discrepancy grows with 


frequency suggest that the main source of error here is 
that caused by the PN truncation. 

We can analyze this further by looking at how the am¬ 
plitude ratio varies with the PN order used to compute 
the modes. Fig.[8]shows the ratio ri m for the (2, 2) mode 
and the (3, 3) mode (which exhibited different behaviors 
in Fig. [7]) computed between the SXS extrapolated NR 
modes and the PN ones including different PN correc¬ 
tions. We recall here that the (2, 2) mode amplitude is 
known up to 3.5 PN while all the other modes are known 
up to 3PN. [We note however that while this paper was in 
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FIG. 7. Ratio of the NR/PN amplitudes as a function of matching frequency Muo for a q = 8 non-spinning binary for several 
NR waveforms matched to Tl. The simulations correspond to r = (100,133,190,266, 307)M and extrapolated SXS data (in 
color and downwards) and BAM r = (60, 80,100 )M (black and upwards). 
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FIG. 8. Ratio of NR and PN amplitudes of q = 8 non-spinning SXS NR data extrapolated at IV = 4 matched to Tl including 
different PN amplitude corrections as a function of Mu q. We show the (2,2) (left) and (3,3) (right) modes. 


preparation, Faye et al. m computed the 3.5PN (non¬ 
spinning) contribution to the (3, 3) mode. This correction 
has not been taken into account in the rest of this paper 
(it would only marginally affect the phasing results shown 
in the left panel of Fig. |5] and the right panels of Figures 
mm and [IT] as well as the match plots in Sec. |VIII[ 
the qualitative behavior remaining unchanged), and in 


particular not in Fig. [7] but we do add it in Fig. [8]where 
it does have an important effect since it cures the large 
discrepancy observed in Fig. [7] and brings the agreement 
between PN and NR to the level of ~ 2% (the ampli¬ 
tude plots Figs. [5] (right) and [b] (bottom left) would also 
noticeably improve if this correction was added). This 
significantly better agreement after improving the PN de- 
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scription of the mode provides additional evidence that 
PN truncation was the dominant source of error here.] 
As usual, studying the PN truncation error is more dif¬ 
ficult than studying for instance the convergence with 
extraction radius because we do not know a priori how 
the series converges and therefore we cannot extrapolate 
it. In particular, the convergence is not necessarily mono¬ 
tonic and consecutive corrections can have very different 
magnitudes (in part because they originate from or mix 
several physical effects) so the fact that adding a given 
correction barely changes the result gives no guarantee 
that the next one will also be negligible. In other words, 
one cannot estimate the truncation error by comparing 
the result at n PN and the one at (n+1/2) PN as trivially 
illustrated by the fact that the 2PN and 2.5PN curves of 
the right panel of Fig. [8] are superimposed but differ sig¬ 
nificantly from the 3PN curve. Despite these caveats, the 
spread between several consecutive curves gives a sense of 
how much PN truncation affects the result and the com¬ 
parison between the left and right panels suggests that 
this error is significantly larger for the (3,3) mode than 
for the (2,2) mode where moreover the highest PN orders 
agree well with NR. We observed a similar spread for the 
(4,3), (4,4) and (5,5) modes while the (2,1) and (3,2) 
modes exhibited a behavior similar to the (2,2) mode. 
Note that while all higher modes are known up to ab¬ 
solute 3PN order (with the convention that the leading 
order of the (2,2) amplitude is Newtonian), their lead¬ 
ing order is {£ — 2)/2 PN for even values of £ + m and 
(£ — 1) /2 PN for odd values of i + m so the number of rel¬ 
ative corrections actually known for each mode of course 
varies. However, we highlight that the differences that 
we observe are not a consequence of a relative higher or¬ 
der knowledge for some modes: both the (2,1) and the 
(3, 3) modes for instance are used with 2.5 PN relative 
precision and show a very different behavior. The lack of 
apparent error systematics regarding PN orders or modes 
is not untypical for PN results, where contributions at 
different orders often come from very different physical 
effects, thus their magnitude is hard to anticipate. 


B. Comparison of PN and NR phases 


We now move to identifying the origin of the residual 
phasing errors defined in Eq. (27). As discussed in 
Sec.[VTJ these quantify the discrepancies between PN and 
NR in the phase difference between the (£, m) mode and 
the (2, 2) mode, i.e. clephasing errors in addition to the 
difference in tracking the orbital phase of the system. 
Indeed, we can rewrite Eq. (271 using Eq. (25) to obtain 


= (« - J «*) - (ffi - “«*) (30) 

where the phases have to be taken at the matching point 
corresponding to the hybridization frequency wo- Note 
that in principle the rhs of Eq. (301 should also contain 


an additional term (2 n'm + (2 — m)n)7 t/2 ) originating 


from convention differences which we assume for simplic¬ 
ity have been reabsorbed in the definition of one of the 
two waveforms. In other words, we adjust the conven¬ 
tions (but only by integer factors of 7t/2 for tpo and of 7r 
for ipo) of say the PN waveform so that the e^ m vanish 
in the limit where both the NR and the PN waveform 
would be infinitely accurate. 

The ee :m are more intricate quantities to study than 
the r£ m since they not only involve PN and NR but also 
two different modes. As a first step, we find it useful to 
focus on the simpler quantities 

A f m (t) = <™W-y^ a W, (31) 


where X is NR or PN, which only involve either the NR 
or PN waveform. Note that A f m is insensitive to a redef¬ 
inition of the angle <p, i.e. with our previous notation, to 
a change in <po. However, a redefinition ipo —s * ipo + Sip o 
affects (311 as A f m — > A.f m + Sip o(l — m/2). In particu¬ 
lar, we define these quantities with the idea in mind that 
during the inspiral they should vanish as the frequency 
decreases, which implies some particular convention for 
ipo- In the rest of this section, in order to simplify the 
discussion, we assume that all the waveforms adopt this 
convention. Note that the conventions adopted in Arun 
et al m and Blanchet [6] differ from this by a shift of n 


in ip 0 . 

In PN, one finds that the A™ are small but non-zero 
during the inspiral, and vanish in the limit of infinite sep¬ 
aration. Consequently, the PN phase of the (£, m ) mode 
approximately follows the rule ~ m(p m o{t) where 

</>orb is the orbital phase of the system. More precisely, in 
the absence of precession, the deviations to this expres¬ 
sion are due to imaginary coefficients in the mode ampli¬ 
tudes which only appear at high PN orders for the modes 
we consider (see Eq. (327) of [B] for the non-spinning 
mode amplitudes; note that the spin corrections to the 
mode amplitudes that are currently known, i.e. up to 
2PN, contain no such complex correction). 

We find excellent agreement between extrapolated NR 
and PN regarding the Af jT71 , i.e. in our NR data we find 
consistent small but non-zero values of A^. However, 
errors resulting from finite extraction radius can be large, 
and in the following we will discuss the dependence of 
the quantities A^ on the extraction radius. Since these 
quantities are functions of time, we first need to align in 
time waveforms extracted at different radii. To facilitate 
the comparison with PN, we pick a reference frequency 
wo early in the waveform and align all the 22 modes with 
say the one with the largest extraction radius using the 
same procedure as described in Sec. 0 Fig - 0 shows the 
(aligned in time) A^ m for the (3,3) and (2,1) modes for 
the q = 8 non-spinning waveform from the SXS catalog 
for different extraction radii as well as for the N = 4 ex¬ 
trapolated numerical waveform and the PN T1 waveform. 
The alignment has been performed at M wq = 0.043, and 
the time coordinate used for the plot has been shifted 
so that t = 0 corresponds to Mw 0 = 0.043 for the ex¬ 
trapolated waveform. For the (2,1) mode, the difference 
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with the extrapolated waveform remains of the order of 
15 — 20 degrees even for the outermost available radii, 
and we observe the same behavior for all the (£, £ — 1) 
modes. On the contrary, for the (3,3) mode (and the 
other (£,£) modes), the larger finite radii curves only dif¬ 
fer from the extrapolated one by of the order of 2 degrees. 
Most importantly, in both cases, the extrapolated wave¬ 
form agrees very well with the PN one, with a typical 
difference of only one degree, illustrating that the main 
source of disagreement in the A^ m comes from the finite 
radius extraction. 

To quantify this further, we studied the asymptotic 
behavior of AA^ m = |A^f — A™ | (averaged in time over 
the interval [1000M, 2000M] to average out oscillations) 
as a function of the extraction radius by fitting it to the 
power law l°(r/ro) n . The best fit values for n and r 0 are 
shown in Table [I] and are consistent with an 1 /r falloff 
rate. 


(Cm) 

(2,1) 

(3,2) 

(3,3) 

(4,3) 

(4,4) 

n 

-0.967 

-1.015 

-0.941 

-1.038 

-0.947 

ro 

3199 

4215 

293 

4182 

598 


TABLE I. Results for (ro,n) for the fits AA £ m = l°(r/ro) n 
for the case q = 8 non-spinning SXS NR data matched to Tl. 
The values suggest an asymptotic 1/r-falloff in the inspiral 
region, where we hybridize. 


We show the values of the e^ m as a function of fre¬ 
quency in Fig. [TO] As expected from our analysis of the 
Ag m , the errors are again dominated by the finite wave 
extraction radius. The largest radii still differ signifi¬ 
cantly from the extrapolated waveforms, for which the 
agreement between NR and PN is of the order of only 
1°. Again, the behavior of the (2,1) mode is typical for 
the (£,£—1) modes, while the (£,£) modes behave similar 
to the (3,3) shown here. We note that the finite radius 
errors are strongly dependent on the gauge conditions 
used in the numerical relativity code. 

We have also analysed the non-spinning q = (3,4,6) 
and the aligned-spin (q = 3, x = ±0.5) cases. We find 
that for systems with different mass ratios and spins 
the quantities (r^ m , A^ m , ef, m ) behave qualitatively the 
same as function of extraction radius. Furthermore, in 
Fig 11 we see that the values of for all the men¬ 
tioned cases are almost the same for equal extraction ra¬ 
dius along all the studied frequency range. This suggests 
that the influence on the e^ >m of the extraction radius r 
of NR waveforms depends on the specific parameters of 
the simulated systems rather weakly. 


VIII. EFFECT OF NR EXTRACTION RADIUS 
AND EXTRAPOLATION ON THE MATCH 


In the previous sections, we have investigated the dis¬ 
agreement between PN and NR in the region where we 
align and attach them. The discrepancies we observed 
illustrate inaccuracies in both waveforms at the typical 


frequencies where we perform the hybridization which 
will contribute to the global error budget of the hybrid. 
However, the accuracy of the full hybrid is of course 
also affected by the details of how these discrepancies 
are smoothed in the hybridization window and most im¬ 
portantly by the intrinsic error of the PN and the NR 
portions before and after the hybridization region. There, 
the rt m and e^ m do not inform us about the accuracy. For 
instance, it could be that the spurious relative dephas¬ 
ings between the NR modes that we have observed may 
have disappeared around the merger where the higher 
modes are most important for matches. It is therefore 
not possible to directly translate the values we observed 
for the residual phase and amplitude disagreements into 
an overall error of the final hybrid. 

In this section, we take a first step towards quantifying 
the error budget of the full hybrid in terms of quantities 
useful for data analysis applications. As usual, we will 
replace the unanswerable question of how much the hy¬ 
brid differs from the true signal by a study of the typical 
mismatch that one obtains when one varies the differ¬ 
ent ingredients in the procedure. A systematic study, 
which properly addresses the different requirements for 
the detection and parameter estimation problems across 
a significant portion of the black hole parameter space, is 
beyond the scope of this paper. Instead, we restrict the 
scope of this study to understanding the effect of the ex¬ 
traction radius of the NR waves (and their extrapolation) 
on the match. In doing this, we do not claim that this is 
the main source of inaccuracy for the hybrid: using for 
instance a PN approximant that predicts a phase evolu¬ 
tion very different from the NR one (see e.g. T4 in Fig. [4]) 
will certainly lead to larger mismatches. Rather, moti¬ 
vated by the observations of the previous section that 
suggest that extraction radius can have a strong effect 
on the agreement between PN and NR in the late inspi¬ 
ral, we illustrate here how much of an effect it makes on 
the match between full hybrids. Comparing the results 
to other studies in the literature can give a first impres¬ 
sion of the relative contribution of different imprecisions 
for data analysis applications and guide more systematic 
future investigations. 

We illustrate our results with the q = 8 non-spinning 
case from the SXS catalog and use the 2015 early Ad¬ 
vanced LIGO noise curve m- We hybridize with Tl (the 
PN approximant that gives the smallest secular trend be¬ 
tween the PN and the NR (2, 2) mode frequencies for this 
case) hybridized at ojq = 0.073. With this choice, and 
given the length of our blending windows, the NR (2, 2) 
mode covers the entire instrumental band (say starting 
from 20 Hz) for total masses larger than ~ 120M 0 . Here 
we have chosen to hybridize at a relatively large fre¬ 
quency to facilitate comparisons with shorter NR sim¬ 
ulations. We focus for now on the highest resolution 
available (namely Lev. 5) and on the waves extracted at 
r = 100M (innermost radius available) and 307M (out¬ 
ermost one) as well as those extrapolated to null infinity 
using N = 2 and TV = 4 polynomial orders. We denote 
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FIG. 9. Value of for different extraction radii (100,133,190, 266, 307)A/ for SXS q = 8 non-spinning data and PN Tl. We 
show here the 21 and 33 modes. The vertical line denotes the location of the merger. The gray and blue curves correspond to 
PN and extrapolated N = 4 data respectively. 



FIG. 10. ee t rn values for a the q = 8 non-spinning system from the SXS catalog for extraction radii r = (100,133,190, 266, 307 )M 
(in color from top to bottom) and extrapolated N = 2,4 (black, respectively solid and dashed) and Tl. We also show the BAM 
r = 100 case in dashed blue. 




(q,A')=(3,0) 

(q,y)=(4,0) 

(q,^)=(6,0) 

(q,A')=(8,0) 

(q,AT)=(3,+0.5) 

(q,^)=(3,-0.5) 


FIG. 11. er,m values for several SXS catalog systems for r = 100A/ and extrapolated N = 4 data (solid and dashed respectively). 
Non-spinning systems are shown in color while we use black and grey for positive and negative spins respectively. Note that 
the values of ee,m are almost equal for all the systems, specially among the non spinning ones. 
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these waveforms /iioo, ^- 307 , ftjv =2 and /ijv =4 respectively. 

Our hybridization procedure applied to each of these 
numerical waveforms yields a hybrid hx{9, <p, k) (with 
the notation of Eq. (|8l) whose modes are defined in 
Eq. (291 and slightly different values for the shift <p Q 
needed to adjust conventions [53] between PN and NR. 
This is of course a mere consequence of the fact that our 
procedure to infer the differences in conventions is af¬ 
fected by the inaccuracies in the original waveforms (here 
we know that the conventions are the same for all the NR 
waveforms that we consider). Because in the presence of 
higher modes, the dependence on ip is non trivial, one has 
to be careful to compare waveforms at the same physi¬ 
cal sky orientation (or optimize over it depending on the 
application). In the general case, if two hybrids are built 
out of different PN pieces and different NR pieces, com¬ 
paring at the same angle ip makes no sense in principle 
since the definitions of the origin of the azimuthal angle 
are a priori independent. In the present case however, 
ip has the same meaning for all the numerical waveforms 
that we consider (and therefore, with our choice of ap¬ 
plying the ipo rotation to the PN part in Eq. (29), for the 
hybrids) and it makes sense to compute 

maxO [/ii(7r/2,0,0),/i 2 (7r/2,0,0)], (32) 


i.e. the overlap optimized over time-shifts only of the 
two waveforms at the same (source) sky location (9 = 
7 t/2 ,<^ = 0) chosen so that the higher modes contribute 
significantly to the full waveform and at the arbitrary ef¬ 
fective polarization k = 0. This is the quantity plotted 
in plain green in Fig. [12] Note however that the early 
inspiral tails of our hybrids now differ by a shift in ip, i.e. 
despite the fact that we are using the same PN input in 


all our hybrids, the quantity defined in Eq. (321 will not 
go exactly to 1 in the limit of small masses where only 
the PN tail is in band [54], Often in data-analysis appli¬ 
cations, it is also interesting to optimize this match over 
ip, in which case all the subtleties of the exact definition 
of the azimuthal angle of the hybrid become irrelevant. 
We therefore also plot in dashed-green the quantity 

maxO [/ii( 77 / 2 ,0,0), /i2(tt/2, ¥>, 0)] (33) 

to,<p 

in the figures below. 

Fig. 12 (beware of the different scale for the top left 
panel) displays the result of this study for various cou¬ 
ples of waveforms. In each panel, we additionally show 
the overlap between the individual modes (in black and 
red) optimized over time and phase to check to what 
extent the modes of both hybrids agree if we allow our¬ 
selves to align them one by one independently. We first 
focus on the top panel. In order to interpret these plots, 
one should keep in mind that in the large mass limit, 
the hybridization window is pushed to lower frequen¬ 
cies than those accessible to the instrument and this be¬ 
comes a pure NR comparison. In this region, the typ¬ 
ical mismatch between hx —4 and /iioo is of a few 0 . 1 % 


and goes below 0.1% for / 1307 . As we move to smaller 
masses and the hybridization region enters the instru¬ 
mental band, the match degrades and reaches a mini¬ 
mum around 40Mg: while the mismatch there is above 
1% in the R = 100 vs N = 4 case (this gets reduced by 
a factor ~ 2 after (^-optimization), it remains as low as 
0.2% when using R = 307. At even lower masses, the 
comparison becomes dominated by the PN tail which is 
identical for both waveforms up to the ip shift discussed 
above and the match grows again (the optimized one go¬ 
ing to 1 exactly in the low mass limit). 

From this comparison where the N = 4 waveforms 
has been used as a reference, one can argue that for the 
NR data set studied here, extracting the waves in the 
SpEC code at radii of 0(300 M) (as typically available 
in the SXS catalog) is sufficient to control the error to 
the ~ 0 . 1 % level in terms of mismatch. 


To check the effect of extrapolation orders on this 
study, in the lower right panel we reproduce the upper 
right one but use N = 2 instead of N = 4. While the gen¬ 
eral behavior and scale remains essentially unchanged, we 
note that the (4,4) and (4,3) modes which were signif¬ 
icantly disagreeing at high frequencies and dominating 
the total mismatch between N = 4 and R = 307 now 
give much higher matches at high masses (note also that 
after optimization over ip the matches between the hy¬ 
brid become very close to 1). This is consistent with 
the Ref. [21], where it is also observed that different or¬ 
ders may show varying performance during the evolution, 
and higher orders may in particular be problematic dur¬ 
ing merger-ringdown. Except for these discrepancies at 
high mass which come from the presence of unphysical 
features in these two N = 4 modes and which remain on 
the order of ~ 0.1%, the order of extrapolation therefore 
does not seem to affect our previous conclusion, as also 
illustrated by the lower left panel in which we directly 
compare the hybrids built with N = 4 and N = 2. 

Finally, we perform a similar study to compare the 
effect of numerical resolution on the match to the ef¬ 
fect from different extraction radii. In Fig. 13 we plot 
matches as in Fig. |12[ but using hybrids constructed from 
numerical waveforms produced at different numerical res¬ 
olutions instead of extraction radii. For the faster than 
polynomial convergence exhibited by spectral codes such 
as SpEC one typically quotes the difference between the 
highest and next-highest resolution as an error estimate, 
see e.g. [Sj for a comparative discussion of the error anal¬ 
ysis performed in spectral and finite-difference numerical 
relativity codes. Here we conservatively use the highest 
available resolution level (5) and a resolution two levels 
coarser (3), and we find mismatches smaller or at com¬ 
parable level that the mismatches shown in Fig. [12] re¬ 
sulting from finite extraction radius. In both cases, the 
mismatches are certainly not larger than what can be ex¬ 
pected in waveform models such as HU EU EUI22- in 
addition, any actual signal searches or parameter estima¬ 
tion calculations will be based on discrete or continuous 
waveform families and matches will effectively be opti- 
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N=4 vs. r=100M N=4vs.r=307M 




M[M 0 ] 

N=4 vs. N=2 



M[M 0 ] 

N=2 vs. r=307M 



M[M 0 ] 


M[M 0 ] 


FIG. 12. Match of individual modes (optimized over phase and time) and full waveforms at (9, ip, i/j) = (7t/ 2,0,0) optimized 
only over time (in solid green) and time and phase (dashed green) for q = 8 non-spinning hybrid waveforms built out of T1 
and SXS NR waveforms produced during a single simulation but either extracted at finite radius r = 100 M and r = 307 M or 
extrapolated with polynomial orders N = 2 and N = 4. 


mized over all physical parameters. 


r=307M: res=5 vs. res=3 



FIG. 13. Match as in Fig. |12| but between hybrids con¬ 
structed from numerical waveforms produced at different nu¬ 
merical resolutions instead of extraction radii. We use the 
highest available resolution level (5) and a lower resolution 
level (3) for R = 307 M. 


IX. DISCUSSION 

Gravitational wave searches based on the data from 
ground based detectors have to date been based on tem¬ 
plate banks which only contain the dominant (l, |m|) = 
(2, 2) spherical harmonic modes of the signal. Several 
authors have recently started to evaluate the impact of 
this restriction on detection and parameter estimation 
[HEMS]. Such studies have used multi-mode effective- 
one-body models and hybrid waveforms directly con¬ 
structed from the numerical data and an inspiral model. 
In addition, such multi-mode hybrids will facilitate the 
extension of phenomenological musi] and similar wave¬ 
form models to the multi-mode case. 

The goal of the present paper has been to achieve a 
better understanding of the construction and properties 
of such hybrids for complete waveforms. To this end, 
we have first compared the subdominant modes of non- 
precessing (and mostly non-spinning) NR data sets with 
PN expressions. We have first described the ambiguities 
in such a comparison, and described how to determine 
the three parameters that parameterise different conven¬ 
tions in the literature and available data sets: a time 
shift r, rotation around the orbital angular momentum 
by an angle ipo, and a polarisation angle tpo- For non- 
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precessing systems it can be assumed that tpo G {0,7r}. 
We have given a prescription to determine ipo from nu¬ 
merical data, although in principle it should be deducible 
from the description of a NR code or PN calculation. 

The presence of PN truncation errors and of numerical 
errors in NR waveform descriptions results in amplitude 
and phase discrepancies during the late inspiral (where 
we can compare them), which cannot be compensated 
by the choice of t, ipo and possibly ipo, and thus some 
choices have to be made to determine the values of these 
quantities for a given multi-mode waveform. In our con¬ 
struction we first determine the time shift r from the 
(2,2) mode alone, and then determine the value of ipo 
from two modes, one of which is typically (2,2). We have 
then studied the deviations between the PN result and 
NR data sets as a function of matching time. 

Regarding NR amplitude errors, we find a strong de¬ 
pendence on the numerical code used, which we attribute 
to the difference in coordinate gauge conditions. With 
the exception of the (3,2)-mode, e.g. r = 100 BAM data 
are significantly closer to the extrapolated result than 
the corresponding SXS r = 100 curve. We also note 
that, for any given finite extraction radius, the error be¬ 
comes larger at lower frequencies. This is the expected 
consequence of the fact that as the frequency increases 
(or equivalently as their wavelength decreases), the wave 
zone (defined by r > A) extends to smaller radii. We 
find good agreement between post-Newtonian amplitude 
and those extrapolated to infinity from several extraction 
radii in a numerical relativity calculation for the (2,1), 
(2, 2) and (3, 2) modes, but significantly larger errors of 
up to a few 10% for the modes (3,3), (4,3), (4,4) and 
(5,5), and we find that these deviations are consistent 
with the spread between PN results at different orders. 

Regarding the phase differences Q im , we find excellent 
agreement on the order of 1° of PN results with numerical 
relativity waveforms extrapolated to null infinity, while 
deviations at finite radius can be as large as tens of de¬ 
grees for (£, £ — 1) modes. On the contrary, for the (£, £) 
modes, the larger finite radii curves only differ from the 
extrapolated one by of the order of 2 degrees. In par¬ 
ticular, we find that the influence on the e^ >m values of 
the extraction radius r of NR waveforms depends on the 
specific parameters of the simulated systems only rather 
weakly. Our results imply that during the inspiral the 
complex part of the standard PN waveform amplitudes 
can be taken as the correct value for practical purposes. 


This yields in particular a convenient test for finite radius 
errors in numerical relativity (since a practically exact 
result is known), and may serve to determine favorable 
coordinate gauge conditions for wave extraction. 

A systematic study of the effect of the phase and ampli¬ 
tude errors which we discuss on gravitational wave data 
analysis, both regarding the detection problem and pa¬ 
rameter estimation is beyond the scope of this paper. In¬ 
stead we have performed simple match calculations with 
the initial and design sensitivity advanced LIGO noise 
curve, evaluating the match between waveforms result¬ 
ing from different extraction radii or numerical resolu¬ 
tions. This simplistic study can serve for comparisons 
due to mismatches resulting from other effects, such as 
the choice of PN approximant, or waveform modelling 
errors, and thus give a first impression of the relative 
importance of different types of waveform imprecisions, 
and help guide more detailed investigations. This part 
of our study overlaps with other investigations such as 
ED Hi El] regarding the errors in finite radius and ex¬ 
trapolated waveforms. Our results appear consistent 
with previous work, but add the aspect of considering 
full hybrids and investigating how the match varies when 
the hybridization region is in band and puts the focus on 
errors in higher modes. 

In Sec. |VII| we have seen in particular that the ex¬ 
traction radius can lead to significant phase and ampli¬ 
tude errors. Corresponding mismatches as discussed in 
Sec. VIII are however roughly at the level of 0.1% for 
the cases we have considered. Mismatches at this level 
appear at least negligible for GW searches. We con¬ 
clude with the remark that a systematic understanding 
of multi-mode waveform errors for parameter estimation 
necessarily needs to take into account precession effects. 
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